6. ARIMA Modeling

(6-1) MODEL NON-STATINARY TS

Nao Mimoto - Dept. of Statistics : The University of Akron

TS Class Web Page – R resource page

---

1. Non-stationary time series

Daily Price of Dow Jones

• From Aug. 28 to Dec. 18, 1972 (Brockwell p.153)

See Shumway p145 GNP modeling with ARIMA

  D <- read.csv("https://nmimoto.github.io/datasets/dowj.csv", header=TRUE)
  Dowj  <- ts(D, start=1, freq=1)    # turn D into TS object
  plot(Dowj, type="o")               # plot with both line and circle

• Stationary?

Lake Huron Data

• Annual Level of Lake Huron 1875-1972

  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/lake.csv", header=TRUE)
  Lake  <- ts(D, start=1875,, freq=1)  # turn it into TS object startng year 1875
  plot(Lake, type="o")

• Stationary?

Global Temperature Data

Global mean land-ocean temparature data from Shumway p.4

  library(astsa)     # do install.packages("astsa") if not installed yet
  data(gtemp)        # load data that is inside astsa package
  plot(gtemp, type="o", ylab="Global Temp Deviations")

• Stationary?

---

2. Tests for Stationarity

How do you check if a time seirs is stationary?

• Visual inspection for constant mean, constant variance

• KPSS test (\(H_0\): Stationary)

• ADF test (\(H_0\): Not Stationary)

• PP test (\(H_0\): Not Stationary)

• Fit AR(1), and see if \(\hat \phi_1\) is significantly different from 1.

KPSS test

Kwiatkowski-Phillips-Schmidt-Shin (1992) test for

• Default choice in auto.arima() \[ \left\{ \begin{array}{ll} H_0: X_t \mbox{ is trend stationary}\\ H_A: X_t \mbox{ is not stationary} \\ \end{array} \right. \]

• Large p-value means “Stationary”.

• Decompose the series as \[ X_t \hspace{3mm} = \hspace{3mm} T_t + W_t + Y_t \] where \(T_t\) is deterministic trend, \(W_t\) is random walk, and \(Y_t\) is stationary error.

• Lagrange multiplier test that the random walk has zero variance.

ADF test

Augumented Dickey-Fuller test (Brockwell p. 194)

• “Unit-root” test

• The stationarity hypothesis \[ \left\{ \begin{array}{ll} H_0: Y_t \mbox{ is not stationary} \\ H_A: Y_t \mbox{ is stationary} \end{array} \right. \]

• Is replaced by \[ \left\{ \begin{array}{ll} H_0: Y_t \mbox{ has unit root} \\ H_A: Y_t \mbox{ does not have unit root } \end{array} \right. \]

• Small p-value rejects the null hypothesis of non-stationarity, and means “Stationary”.

• Fit AR(1) to a time series, and test if \(\phi_1\) is significantly different from 1, by \[ \frac{ \hat \phi_1 -1 }{ SE(\hat \phi_1) } \sim N(0,1) \]

PP test

Phillips-Perron test for the null hypothesis that x has a unit root.

• Same Hypothesis as ADF test

• Small p-value means “Stationary”.

• The tests are similar to ADF tests, but they incorporate an automatic correction to the DF procedure to all ow for autocorrelated residuals.

• Calculation of the test statistics is complex.

• The tests usually give the same conclusions as the ADF tests

P-value of (KPSS, ADF, PP) tests

* (LARGE, small, small) \hspace{3mm} $\to$ All three indicating stationarity.

* (small, LARGE, LARGE) \hspace{3mm} $\to$ All three indicating non-stationarity. 

* (anything else) \hspace{3mm} $\to$ Conflicing conclusions, inconclusive. 

• Stationarity.tests() performs all of the three tests.

 #- same file as Randomness.tests()  
 source("https://nmimoto.github.io/477/TS_R-90.txt")  
                         # load Stationarity.tests() from class web page.

---

Ex: Use it on above datasets

  source("https://nmimoto.github.io/477/TS_R-90.txt")  
                          # Load Stationarity.tests() from Course Webpage

  #--- Test it on random sample ---
  X <- rnorm(100)
  plot(X, type="o")

  Stationarity.tests(X)  

## Warning in adf.test(A): p-value smaller than printed p-value

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value greater than printed p-value

##        KPSS  ADF   PP
## p-val:  0.1 0.01 0.01

• Large, small, small \(\to\) stationary. As it should be.

  #--- Lake Huron Data ---
  plot(Lake, type="o") 

  Stationarity.tests(Lake)

## Warning in adf.test(A): p-value smaller than printed p-value

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value greater than printed p-value

##        KPSS  ADF   PP
## p-val:  0.1 0.01 0.01

• small, Large, small \(\to\) conflicting. First two is indicating Non-Stationarity, last one is indicating stationarity.

  #--- Global Temp Data ---
  plot(gtemp, type="o", ylab="Global Temperature Deviations")

  Stationarity.tests(gtemp)

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value smaller than printed p-value

##        KPSS   ADF   PP
## p-val: 0.01 0.706 0.01

• small, Large, Large \(\to\) All three indicating non-stationarity.

---

Ex: Australian Steel Data

  D  <- read.csv("https://nmimoto.github.io/datasets/Steel.csv")
  D1 <- ts(D[,2], start=c(1956,1), freq=12)
  D2 <- window(D1, start=c(1969, 12))

  plot(D1, type="o")

  Stationarity.tests(D1)

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value smaller than printed p-value

##        KPSS   ADF   PP
## p-val: 0.01 0.408 0.01

  plot(D2, type="o")

  Stationarity.tests(D2)

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in pp.test(A): p-value smaller than printed p-value

##        KPSS   ADF   PP
## p-val: 0.01 0.049 0.01

  D4 <- window(D2, end=c(1981,  12))    # Dec 1969 to Dec 1981
  D5 <- window(D2, start=c(1984, 1))    # Feb 1984 to Nov 1993

  Stationarity.tests(D4)

## Warning in adf.test(A): p-value smaller than printed p-value

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value smaller than printed p-value

##        KPSS  ADF   PP
## p-val: 0.01 0.01 0.01

  Stationarity.tests(D5)

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value greater than printed p-value

##        KPSS   ADF   PP
## p-val:  0.1 0.046 0.01

---

3. Modeling Lake HURON data 4 different ways

a) Direct ARMA fit

If we treat it as “stationary” series, we can try to fit ARMA(p,q) series.

  library(forecast)
  source('https://nmimoto.github.io/R/TS-00.txt') 
  
  D <- read.csv("https://nmimoto.github.io/datasets/lake.csv")
  Lake  <- ts(D, start=1875, freq=1)
  
  plot(Lake, type="o")

  layout(matrix(c(1,2), 1, 2))
  acf(Lake)
  pacf(Lake)

  layout(1)

  Stationarity.tests(Lake)

## Warning in kpss.test(A): p-value smaller than printed p-value

##        KPSS  ADF    PP
## p-val: 0.01 0.24 0.025

  # find best ARMA(p,q) by AICc
  Fit1 <- auto.arima(Lake, d=0, stepwise=FALSE, approximation=FALSE)  
              
  Fit1    

## Series: Lake 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    mean
##       0.7665  0.3393  9.1290
## s.e.  0.0773  0.1123  0.3861
## 
## sigma^2 estimated as 0.4784:  log likelihood=-101.09
## AIC=210.18   AICc=210.62   BIC=220.48

  Randomness.tests(Fit1$residuals)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20   JB    SD
## [1,] 0.963 0.952 0.934 0.567 0.641 0.89 0.684

Analysis 1 (direct fit)

• auto.arima() chooses ARMA(1,1) with min AICc.

• AR(1,1) with constant mean was fit directly to data. With observation \(Y_t\),

\[{\large Y_t \hspace{3mm} = \hspace{3mm} \mu + X_t }\] \[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} }\] \[{\large e_t \sim WN(0, \sigma^2) }\]

  Fit1

## Series: Lake 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    mean
##       0.7665  0.3393  9.1290
## s.e.  0.0773  0.1123  0.3861
## 
## sigma^2 estimated as 0.4784:  log likelihood=-101.09
## AIC=210.18   AICc=210.62   BIC=220.48

---

b) Forcing Linear Trend

If we choose that the level is going down linearly, we can try to fit a line with ARMA errors.

  Reg2 <- lm(Lake~time(Lake))
  summary(Reg2)

## 
## Call:
## lm(formula = Lake ~ time(Lake))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.50919 -0.74760 -0.01556  0.75966  2.53409 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 55.278141   7.922614   6.977 4.02e-10 ***
## time(Lake)  -0.024071   0.004119  -5.843 7.16e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.136 on 95 degrees of freedom
## Multiple R-squared:  0.2644, Adjusted R-squared:  0.2566 
## F-statistic: 34.14 on 1 and 95 DF,  p-value: 7.165e-08

  plot(Lake, type="o")                     
  abline(Reg2, col="red")         # overlay regression line in red

  plot(Reg2$residuals, type="o")            # plot of the residuals from regression

  # fit the residuals with ARMA(p,q)
  Fit2 <- auto.arima(Reg2$residuals, d=0, stepwise=FALSE, approximation=FALSE)   
  Fit2

## Series: Reg2$residuals 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##          ar1     ma1
##       0.6671  0.3827
## s.e.  0.0937  0.1135
## 
## sigma^2 estimated as 0.452:  log likelihood=-98.72
## AIC=203.44   AICc=203.7   BIC=211.16

  Randomness.tests(Fit2$resid)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25 ML15  ML20    JB    SD
## [1,] 0.985 0.974 0.969 0.19 0.189 0.782 0.669

Analysis 2 (linear trend)

• Linear trend was fit to \(Y_t\), then the residuals were fitted with ARMA(1,1).

• We observe $ Y_t = a+b t + X_t $ \(X_t\) is ARMA

• In other words,

\[{\large Y_t \hspace{3mm} = \hspace{3mm} a + bt + X_t }\]

\[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} }\]

\[{\large e_t \sim WN(0, \sigma^2) }\]

  summary(Reg2)

## 
## Call:
## lm(formula = Lake ~ time(Lake))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.50919 -0.74760 -0.01556  0.75966  2.53409 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 55.278141   7.922614   6.977 4.02e-10 ***
## time(Lake)  -0.024071   0.004119  -5.843 7.16e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.136 on 95 degrees of freedom
## Multiple R-squared:  0.2644, Adjusted R-squared:  0.2566 
## F-statistic: 34.14 on 1 and 95 DF,  p-value: 7.165e-08

  Fit2

## Series: Reg2$residuals 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##          ar1     ma1
##       0.6671  0.3827
## s.e.  0.0937  0.1135
## 
## sigma^2 estimated as 0.452:  log likelihood=-98.72
## AIC=203.44   AICc=203.7   BIC=211.16

---

c) Box-Cox Differencing (ARIMA modeling)

Let’s look at the difference between the observations: \[ diff(Y) = Y_t - Y_{t-1} \\ \\ \]

  diff.Y <- diff(Lake)      # takes difference between the observations
  plot(diff.Y, type='o')

  layout(matrix(1:2, 1, 2))
  acf(diff.Y)
  pacf(diff.Y)

  layout(1)

  Stationarity.tests(diff.Y)       # See if diff(Y) is stationary  

## Warning in adf.test(A): p-value smaller than printed p-value

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value greater than printed p-value

##        KPSS  ADF   PP
## p-val:  0.1 0.01 0.01

  # Fit the difference of Lake with ARMA 
  Fit3 <- auto.arima(diff.Y, d=0, stepwise=FALSE, approximation=FALSE)    
  Fit3

## Series: diff.Y 
## ARIMA(1,0,2) with zero mean 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.6385  -0.5349  -0.3514
## s.e.  0.1345   0.1445   0.1055
## 
## sigma^2 estimated as 0.4812:  log likelihood=-99.88
## AIC=207.76   AICc=208.2   BIC=218.02

  Randomness.tests(Fit3$residuals)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.981 0.958 0.946 0.551 0.545 0.575 0.681

  #-- Exactly same as above 
  Fit3b <- Arima(Lake, order=c(1,1,2), include.mean=FALSE)    
  Fit3b

## Series: Lake 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.6385  -0.5349  -0.3514
## s.e.  0.1345   0.1445   0.1055
## 
## sigma^2 estimated as 0.4812:  log likelihood=-99.88
## AIC=207.76   AICc=208.2   BIC=218.02

  Randomness.tests(Fit3b$residuals)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.979 0.956 0.942 0.524 0.495 0.577 0.677

Analysis 3 (take difference)

• Difference of \(Y_t\) was taken. The difference seems to be ARMA(1,2).

• In other words, $ Y_t $ is ARIMA(1,1,2), meaning \[{\large \bigtriangledown Y_t = Y_t - Y_{t-1} = X_t }\] \[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} }\] \[{\large e_t \sim WN(0,\sigma^2) }\]

  Fit3b

## Series: Lake 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.6385  -0.5349  -0.3514
## s.e.  0.1345   0.1445   0.1055
## 
## sigma^2 estimated as 0.4812:  log likelihood=-99.88
## AIC=207.76   AICc=208.2   BIC=218.02

---

d) Box-Cox Differencing (ARIMA modeling)

Let’s look at the difference between the observations: \[ diff(Y) = Y_t - Y_{t-1} \\ \\ \]

  diff.Y <- diff(Lake)      # takes difference between the observations
  plot(diff.Y, type='o')

  Randomness.tests(diff(Lake))

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.176 0.175 0.132 0.172 0.117 0.384 0.737

  #-- Exactly same as above 
  Fit4 <- Arima(Lake, order=c(0,1,0), include.mean=FALSE)    
  Fit4

## Series: Lake 
## ARIMA(0,1,0) 
## 
## sigma^2 estimated as 0.5383:  log likelihood=-106.49
## AIC=214.98   AICc=215.02   BIC=217.54

Analysis 4 (take difference)

• Difference of \(Y_t\) was taken. The difference seems to be WN.

• In other words, $ Y_t $ is ARIMA(0,1,0), meaning \[{\large \bigtriangledown Y_t = Y_t - Y_{t-1} = X_t }\] \[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} e_t }\] \[{\large e_t \sim WN(0,\sigma^2) }\]

    Fit4

## Series: Lake 
## ARIMA(0,1,0) 
## 
## sigma^2 estimated as 0.5383:  log likelihood=-106.49
## AIC=214.98   AICc=215.02   BIC=217.54

---

Summary

Three popular test for stationarity is KPSS, ADF, and PP test.

Tests are availabe in R, can be implemented with below command.

  source('https://nmimoto.github.io/R/TS-00.txt')
  Stationarity.tests(X)

p-value of (Large, Small, Small) indicates stationarity.

Classical method tries to identify trend, and tries to removeit by polynomial regression.

Box-Jenkins Method tries to take a difference between today’s data and yesterday’s by applying \(\bigtriangledown=(1-B)\) operator.

Applying B-J method will often times stationarize non-stationary time series.

We can use ARMA(\(p,q\)) model to model the stationarized version of the series.

If you use \(\bigtriangledown\) once, and apply ARMA(1,2), then it is called ARIMA (1,1,2) model.

TS Class Web Page – R resource page